We plot the data and can see that there is no obvious large difference between the versions with high and low debt.
d.both_completed %>%
ggplot(aes(x=time/60, fill=high_debt_version)) +
geom_boxplot() +
labs(
title = "Distribution of time measurements for the different debt levels",
subtitle = "Notice! Log10 x-scale",
x ="Time (min)"
) +
scale_y_continuous(breaks = NULL) +
scale_x_log10() +
scale_fill_manual(
name = "Debt level",
labels = c("High debt", "Low debt"),
values = c("#7070FF", "lightblue"),
guide = guide_legend(reverse = TRUE)
)
d.both_completed %>%
pull(time) %>%
summary()
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 442.0 937.5 1358.0 1842.2 2423.8 7298.0
sprintf("Variance: %.0f", var(pull(d.both_completed, time)))
## [1] "Variance: 1880646"
As the variance is much greater than the mean we will use a negative binomial family that allows us to model the variance separately.
We include high_debt_verison as a predictor in our model as this variable represent the very effect we want to measure. We also include a varying intercept for each individual to prevent the model from learning too much from single participants with extreme measurements.
We iterate over the model until we have sane priors.
time.with <- extendable_model(
base_name = "time",
base_formula = "time ~ 1 + high_debt_version + (1 | session)",
base_priors = c(
prior(normal(0, 1), class = "b"),
prior(normal(7.5, 1), class = "Intercept"),
prior(exponential(1), class = "sd"),
prior(gamma(0.01, 0.01), class = "shape")
),
family = negbinomial(),
data = d.both_completed,
base_control = list(adapt_delta = 0.95)
)
prior_summary(time.with(only_priors= TRUE))
prior_summary(time.with(sample_prior = "only"))
pp_check(time.with(sample_prior = "only"), nsamples = 200) + scale_x_log10()
We check the posterior distribution and can see that the model seems to have been able to fit the data well. Sampling seems to also have worked well as Rhat values are close to 1 and the sampling plots look nice.
pp_check(time.with(), nsamples = 200) + scale_x_log10()
summary(time.with())
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: time ~ 1 + high_debt_version + (1 | session)
## Data: as.data.frame(data) (Number of observations: 44)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~session (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.41 0.15 0.09 0.72 1.00 720 735
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 7.53 0.15 7.24 7.83 1.00 1979
## high_debt_versionfalse -0.18 0.17 -0.50 0.15 1.00 3994
## Tail_ESS
## Intercept 2611
## high_debt_versionfalse 2740
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 3.64 1.08 1.92 6.08 1.00 1213 2177
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(time.with(), ask = FALSE)
# default prior for monotonic predictor
edlvl_prior <- prior(dirichlet(2), class = "simo", coef = "moeducation_level1")
We use loo to check some possible extensions on the model.
loo_result <- loo(
# Benchmark model(s)
time.with(),
# New model(s)
time.with("work_domain"),
time.with("work_experience_programming.s"),
time.with("work_experience_java.s"),
time.with("education_field"),
time.with("mo(education_level)", edlvl_prior),
time.with("workplace_peer_review"),
time.with("workplace_td_tracking"),
time.with("workplace_pair_programming"),
time.with("workplace_coding_standards"),
time.with("scenario"),
time.with("group")
)
loo_result[2]
## $diffs
## elpd_diff se_diff
## time.with("scenario") 0.0 0.0
## time.with("education_field") -1.0 3.8
## time.with("workplace_peer_review") -1.5 2.8
## time.with() -1.5 2.2
## time.with("mo(education_level)", edlvl_prior) -1.7 3.6
## time.with("workplace_pair_programming") -2.0 2.4
## time.with("workplace_coding_standards") -2.3 2.3
## time.with("work_experience_java.s") -2.4 2.2
## time.with("work_experience_programming.s") -2.6 2.3
## time.with("workplace_td_tracking") -2.8 1.9
## time.with("group") -3.5 2.5
## time.with("work_domain") -3.7 4.1
loo_result[1]
## $loos
## $loos$`time.with()`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.6 7.2
## p_loo 14.5 2.8
## looic 735.2 14.3
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 26 59.1% 886
## (0.5, 0.7] (ok) 14 31.8% 331
## (0.7, 1] (bad) 4 9.1% 53
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("work_domain")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -369.7 6.4
## p_loo 16.6 2.5
## looic 739.5 12.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 27 61.4% 656
## (0.5, 0.7] (ok) 15 34.1% 135
## (0.7, 1] (bad) 1 2.3% 30
## (1, Inf) (very bad) 1 2.3% 33
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("work_experience_programming.s")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -368.7 7.3
## p_loo 15.6 3.1
## looic 737.3 14.5
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 27 61.4% 799
## (0.5, 0.7] (ok) 12 27.3% 258
## (0.7, 1] (bad) 5 11.4% 25
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("work_experience_java.s")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -368.4 7.1
## p_loo 15.0 2.8
## looic 736.9 14.1
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 30 68.2% 425
## (0.5, 0.7] (ok) 12 27.3% 388
## (0.7, 1] (bad) 2 4.5% 45
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("education_field")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.1 6.2
## p_loo 12.9 2.0
## looic 734.1 12.3
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 28 63.6% 749
## (0.5, 0.7] (ok) 13 29.5% 286
## (0.7, 1] (bad) 3 6.8% 106
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("mo(education_level)", edlvl_prior)`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.8 6.9
## p_loo 13.2 2.6
## looic 735.6 13.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 26 59.1% 685
## (0.5, 0.7] (ok) 13 29.5% 672
## (0.7, 1] (bad) 4 9.1% 83
## (1, Inf) (very bad) 1 2.3% 19
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("workplace_peer_review")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.6 6.4
## p_loo 12.8 2.3
## looic 735.2 12.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 28 63.6% 557
## (0.5, 0.7] (ok) 15 34.1% 326
## (0.7, 1] (bad) 1 2.3% 47
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("workplace_td_tracking")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -368.9 8.0
## p_loo 16.1 3.8
## looic 737.8 16.0
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 22 50.0% 895
## (0.5, 0.7] (ok) 18 40.9% 405
## (0.7, 1] (bad) 3 6.8% 167
## (1, Inf) (very bad) 1 2.3% 5
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("workplace_pair_programming")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -368.0 7.1
## p_loo 15.0 2.9
## looic 736.1 14.3
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 27 61.4% 805
## (0.5, 0.7] (ok) 13 29.5% 320
## (0.7, 1] (bad) 3 6.8% 43
## (1, Inf) (very bad) 1 2.3% 34
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("workplace_coding_standards")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -368.3 7.1
## p_loo 14.8 2.9
## looic 736.7 14.1
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 30 68.2% 747
## (0.5, 0.7] (ok) 7 15.9% 600
## (0.7, 1] (bad) 6 13.6% 60
## (1, Inf) (very bad) 1 2.3% 14
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("scenario")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -366.1 8.0
## p_loo 17.3 3.7
## looic 732.1 16.1
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 26 59.1% 786
## (0.5, 0.7] (ok) 12 27.3% 328
## (0.7, 1] (bad) 5 11.4% 51
## (1, Inf) (very bad) 1 2.3% 14
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("group")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -369.5 6.8
## p_loo 16.1 2.8
## looic 739.0 13.6
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 23 52.3% 539
## (0.5, 0.7] (ok) 18 40.9% 343
## (0.7, 1] (bad) 3 6.8% 34
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
loo_result <- loo(
# Benchmark model(s)
time.with(),
time.with("scenario"),
time.with("education_field"),
time.with("workplace_peer_review"),
time.with("mo(education_level)", edlvl_prior),
# New model(s)
time.with(c("scenario", "education_field")),
time.with(c("scenario", "mo(education_level)"), edlvl_prior),
time.with(c("scenario", "workplace_peer_review")),
time.with(c("education_field", "mo(education_level)"), edlvl_prior),
time.with(c("education_field", "workplace_peer_review")),
time.with(c("mo(education_level)", "workplace_peer_review"), edlvl_prior)
)
loo_result[2]
## $diffs
## elpd_diff
## time.with(c("scenario", "mo(education_level)"), edlvl_prior) 0.0
## time.with(c("scenario", "education_field")) 0.0
## time.with("scenario") -0.3
## time.with(c("scenario", "workplace_peer_review")) -0.8
## time.with("education_field") -1.3
## time.with("workplace_peer_review") -1.8
## time.with() -1.8
## time.with(c("education_field", "workplace_peer_review")) -1.9
## time.with("mo(education_level)", edlvl_prior) -2.0
## time.with(c("mo(education_level)", "workplace_peer_review"), edlvl_prior) -2.5
## time.with(c("education_field", "mo(education_level)"), edlvl_prior) -2.5
## se_diff
## time.with(c("scenario", "mo(education_level)"), edlvl_prior) 0.0
## time.with(c("scenario", "education_field")) 2.1
## time.with("scenario") 2.5
## time.with(c("scenario", "workplace_peer_review")) 2.0
## time.with("education_field") 3.3
## time.with("workplace_peer_review") 2.7
## time.with() 3.0
## time.with(c("education_field", "workplace_peer_review")) 3.1
## time.with("mo(education_level)", edlvl_prior) 2.4
## time.with(c("mo(education_level)", "workplace_peer_review"), edlvl_prior) 2.5
## time.with(c("education_field", "mo(education_level)"), edlvl_prior) 2.9
loo_result[1]
## $loos
## $loos$`time.with()`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.6 7.2
## p_loo 14.5 2.8
## looic 735.2 14.3
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 26 59.1% 886
## (0.5, 0.7] (ok) 14 31.8% 331
## (0.7, 1] (bad) 4 9.1% 53
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("scenario")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -366.1 8.0
## p_loo 17.3 3.7
## looic 732.1 16.1
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 26 59.1% 786
## (0.5, 0.7] (ok) 12 27.3% 328
## (0.7, 1] (bad) 5 11.4% 51
## (1, Inf) (very bad) 1 2.3% 14
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("education_field")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.1 6.2
## p_loo 12.9 2.0
## looic 734.1 12.3
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 28 63.6% 749
## (0.5, 0.7] (ok) 13 29.5% 286
## (0.7, 1] (bad) 3 6.8% 106
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("workplace_peer_review")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.6 6.4
## p_loo 12.8 2.3
## looic 735.2 12.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 28 63.6% 557
## (0.5, 0.7] (ok) 15 34.1% 326
## (0.7, 1] (bad) 1 2.3% 47
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("mo(education_level)", edlvl_prior)`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.8 6.9
## p_loo 13.2 2.6
## looic 735.6 13.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 26 59.1% 685
## (0.5, 0.7] (ok) 13 29.5% 672
## (0.7, 1] (bad) 4 9.1% 83
## (1, Inf) (very bad) 1 2.3% 19
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("scenario", "education_field"))`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -365.8 6.8
## p_loo 15.6 2.5
## looic 731.7 13.5
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 25 56.8% 435
## (0.5, 0.7] (ok) 14 31.8% 425
## (0.7, 1] (bad) 5 11.4% 50
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("scenario", "mo(education_level)"), edlvl_prior)`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -365.8 7.3
## p_loo 15.2 2.9
## looic 731.6 14.6
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 28 63.6% 391
## (0.5, 0.7] (ok) 13 29.5% 129
## (0.7, 1] (bad) 3 6.8% 49
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("scenario", "workplace_peer_review"))`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -366.6 7.5
## p_loo 16.6 3.3
## looic 733.1 14.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 29 65.9% 457
## (0.5, 0.7] (ok) 10 22.7% 172
## (0.7, 1] (bad) 4 9.1% 177
## (1, Inf) (very bad) 1 2.3% 23
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("education_field", "mo(education_level)"), edlvl_prior)`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -368.3 6.5
## p_loo 14.1 2.3
## looic 736.7 12.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 22 50.0% 628
## (0.5, 0.7] (ok) 17 38.6% 629
## (0.7, 1] (bad) 4 9.1% 62
## (1, Inf) (very bad) 1 2.3% 36
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("education_field", "workplace_peer_review"))`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.7 5.8
## p_loo 12.7 1.8
## looic 735.4 11.6
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 27 61.4% 512
## (0.5, 0.7] (ok) 12 27.3% 501
## (0.7, 1] (bad) 5 11.4% 94
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("mo(education_level)", "workplace_peer_review"), edlvl_prior)`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -368.3 6.8
## p_loo 13.5 2.6
## looic 736.5 13.5
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 24 54.5% 525
## (0.5, 0.7] (ok) 16 36.4% 521
## (0.7, 1] (bad) 3 6.8% 66
## (1, Inf) (very bad) 1 2.3% 16
## See help('pareto-k-diagnostic') for details.
loo_result <- loo(
# Benchmark model(s)
time.with(),
time.with("scenario"),
time.with(c("scenario", "education_field")),
time.with(c("scenario", "mo(education_level)"), edlvl_prior),
# New model(s)
time.with(c("scenario", "mo(education_level)", "education_field"), edlvl_prior)
)
loo_result[2]
## $diffs
## elpd_diff
## time.with(c("scenario", "mo(education_level)"), edlvl_prior) 0.0
## time.with(c("scenario", "education_field")) 0.0
## time.with("scenario") -0.3
## time.with(c("scenario", "mo(education_level)", "education_field"), edlvl_prior) -0.7
## time.with() -1.8
## se_diff
## time.with(c("scenario", "mo(education_level)"), edlvl_prior) 0.0
## time.with(c("scenario", "education_field")) 2.1
## time.with("scenario") 2.5
## time.with(c("scenario", "mo(education_level)", "education_field"), edlvl_prior) 1.4
## time.with() 3.0
loo_result[1]
## $loos
## $loos$`time.with()`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -367.6 7.2
## p_loo 14.5 2.8
## looic 735.2 14.3
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 26 59.1% 886
## (0.5, 0.7] (ok) 14 31.8% 331
## (0.7, 1] (bad) 4 9.1% 53
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with("scenario")`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -366.1 8.0
## p_loo 17.3 3.7
## looic 732.1 16.1
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 26 59.1% 786
## (0.5, 0.7] (ok) 12 27.3% 328
## (0.7, 1] (bad) 5 11.4% 51
## (1, Inf) (very bad) 1 2.3% 14
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("scenario", "education_field"))`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -365.8 6.8
## p_loo 15.6 2.5
## looic 731.7 13.5
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 25 56.8% 435
## (0.5, 0.7] (ok) 14 31.8% 425
## (0.7, 1] (bad) 5 11.4% 50
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("scenario", "mo(education_level)"), edlvl_prior)`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -365.8 7.3
## p_loo 15.2 2.9
## looic 731.6 14.6
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 28 63.6% 391
## (0.5, 0.7] (ok) 13 29.5% 129
## (0.7, 1] (bad) 3 6.8% 49
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## $loos$`time.with(c("scenario", "mo(education_level)", "education_field"), edlvl_prior)`
##
## Computed from 4000 by 44 log-likelihood matrix
##
## Estimate SE
## elpd_loo -366.5 6.9
## p_loo 16.1 2.6
## looic 733.0 13.7
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 29 65.9% 307
## (0.5, 0.7] (ok) 11 25.0% 109
## (0.7, 1] (bad) 3 6.8% 115
## (1, Inf) (very bad) 1 2.3% 12
## See help('pareto-k-diagnostic') for details.
We pick some of our top performing models as candidates and inspect them closer.
The candidate models are named and listed in order of complexity.
We select the simplest model as a baseline.
time0 <- brm(
"time ~ 1 + high_debt_version + (1 | session)",
prior = c(
prior(normal(0, 1), class = "b"),
prior(normal(7.5, 1), class = "Intercept"),
prior(exponential(1), class = "sd"),
prior(gamma(0.01, 0.01), class = "shape")
),
family = negbinomial(),
data = as.data.frame(d.both_completed),
control = list(adapt_delta = 0.95),
file = "fits/time0",
file_refit = "on_change",
seed = 20210421
)
summary(time0)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: time ~ 1 + high_debt_version + (1 | session)
## Data: as.data.frame(d.both_completed) (Number of observations: 44)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~session (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.42 0.15 0.09 0.71 1.00 784 991
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 7.54 0.15 7.25 7.85 1.00 2527
## high_debt_versionfalse -0.18 0.17 -0.51 0.17 1.00 4094
## Tail_ESS
## Intercept 2316
## high_debt_versionfalse 2711
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 3.67 1.10 1.89 6.13 1.00 1114 1722
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
ranef(time0)
## $session
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## 6033d69a5af2c702367b3a95 -0.11275159 0.2860867 -0.685643550 0.4582609
## 6033d90a5af2c702367b3a96 0.30118528 0.2841025 -0.192779725 0.9023268
## 6034fc165af2c702367b3a98 0.31146058 0.2897695 -0.202207675 0.9059238
## 603500725af2c702367b3a99 -0.50158676 0.3736261 -1.239848500 0.1540829
## 603f97625af2c702367b3a9d -0.17573161 0.2929061 -0.763236875 0.3805710
## 603fd5d95af2c702367b3a9e 0.23260419 0.2826572 -0.282267400 0.8325673
## 60409b7b5af2c702367b3a9f 0.36297895 0.2969008 -0.149718000 0.9718112
## 604b82b5a7718fbed181b336 -0.32887367 0.3349227 -1.005135250 0.2731784
## 6050c1bf856f36729d2e5218 0.11487089 0.2748997 -0.388343300 0.6867415
## 6050e1e7856f36729d2e5219 0.07194507 0.2762681 -0.465090150 0.6481566
## 6055fdc6856f36729d2e521b -0.07517223 0.2752623 -0.622037225 0.4849628
## 60589862856f36729d2e521f 0.04470959 0.2798324 -0.476523000 0.6353896
## 605afa3a856f36729d2e5222 -0.07470695 0.2916273 -0.648323775 0.5008381
## 605c8bc6856f36729d2e5223 -0.35466018 0.3355242 -1.037378000 0.2627359
## 605f3f2d856f36729d2e5224 -0.11281797 0.2921469 -0.697781400 0.4773815
## 605f46c3856f36729d2e5225 -0.24574637 0.3096026 -0.862547675 0.3286244
## 60605337856f36729d2e5226 0.22756347 0.2824954 -0.265445925 0.8269080
## 60609ae6856f36729d2e5228 0.31306097 0.2932300 -0.203570125 0.9317685
## 6061ce91856f36729d2e522e -0.04088890 0.2757911 -0.587677450 0.5242268
## 6061f106856f36729d2e5231 -0.36697510 0.3340662 -1.055116250 0.2474501
## 6068ea9f856f36729d2e523e 0.65940365 0.3449495 0.003345279 1.3392622
## 6075ab05856f36729d2e5247 -0.30062645 0.3269950 -0.971954650 0.3085037
plot(time0, ask = FALSE)
pp_check(time0, nsamples = 200) + scale_x_log10()
We select the best performing model with one variable.
time1 <- brm(
"time ~ 1 + high_debt_version + scenario + (1 | session)",
prior = c(
prior(normal(0, 1), class = "b"),
prior(normal(7.5, 1), class = "Intercept"),
prior(exponential(1), class = "sd"),
prior(gamma(0.01, 0.01), class = "shape")
),
family = negbinomial(),
data = as.data.frame(d.both_completed),
control = list(adapt_delta = 0.95),
file = "fits/time1",
file_refit = "on_change",
seed = 20210421
)
summary(time1)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: time ~ 1 + high_debt_version + scenario + (1 | session)
## Data: as.data.frame(d.both_completed) (Number of observations: 44)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~session (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.46 0.14 0.17 0.73 1.00 972 1288
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 7.65 0.16 7.34 7.97 1.00 2512
## high_debt_versionfalse -0.15 0.16 -0.46 0.19 1.00 4702
## scenariotickets -0.30 0.16 -0.62 0.03 1.00 4391
## Tail_ESS
## Intercept 2862
## high_debt_versionfalse 2577
## scenariotickets 3089
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 4.24 1.29 2.17 7.17 1.00 1212 2108
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
ranef(time1)
## $session
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## 6033d69a5af2c702367b3a95 -0.131157520 0.2894703 -0.69790085 0.4472849
## 6033d90a5af2c702367b3a96 0.345975163 0.2948088 -0.17240543 0.9621481
## 6034fc165af2c702367b3a98 0.312615721 0.2890198 -0.21111305 0.9314271
## 603500725af2c702367b3a99 -0.580967048 0.3753131 -1.29102175 0.1374856
## 603f97625af2c702367b3a9d -0.178034878 0.2968035 -0.75711637 0.3855893
## 603fd5d95af2c702367b3a9e 0.216871452 0.2843529 -0.30592422 0.8045455
## 60409b7b5af2c702367b3a9f 0.434482730 0.3044064 -0.10936117 1.0527390
## 604b82b5a7718fbed181b336 -0.384507154 0.3158815 -0.99940722 0.2023649
## 6050c1bf856f36729d2e5218 0.178131596 0.2899054 -0.37398580 0.7811601
## 6050e1e7856f36729d2e5219 0.084611385 0.2835651 -0.44492085 0.6639378
## 6055fdc6856f36729d2e521b -0.067260531 0.2879363 -0.62157198 0.5123126
## 60589862856f36729d2e521f 0.007323255 0.2802704 -0.53228525 0.5927170
## 605afa3a856f36729d2e5222 -0.090608450 0.2803041 -0.63434492 0.4693394
## 605c8bc6856f36729d2e5223 -0.424604250 0.3286785 -1.06578100 0.2004248
## 605f3f2d856f36729d2e5224 -0.084178630 0.2924085 -0.63096325 0.5377237
## 605f46c3856f36729d2e5225 -0.274337832 0.2973105 -0.84274730 0.2896524
## 60605337856f36729d2e5226 0.263460224 0.2959967 -0.27361325 0.8840842
## 60609ae6856f36729d2e5228 0.362073268 0.2928321 -0.17288075 0.9783804
## 6061ce91856f36729d2e522e -0.054703919 0.2956666 -0.63568705 0.5367869
## 6061f106856f36729d2e5231 -0.390416491 0.3256629 -1.03886450 0.2198758
## 6068ea9f856f36729d2e523e 0.804105705 0.3408505 0.07752689 1.4607338
## 6075ab05856f36729d2e5247 -0.357241946 0.3189408 -0.98133657 0.2437772
plot(time1, ask = FALSE)
pp_check(time1, nsamples = 200) + scale_x_log10()
We select the best performing model with two variables.
time2 <- brm(
"time ~ 1 + high_debt_version + scenario + mo(education_level) + (1 | session)",
prior = c(
prior(normal(0, 1), class = "b"),
prior(normal(7.5, 1), class = "Intercept"),
prior(exponential(1), class = "sd"),
prior(gamma(0.01, 0.01), class = "shape"),
prior(dirichlet(2), class = "simo", coef = "moeducation_level1")
),
family = negbinomial(),
data = as.data.frame(d.both_completed),
control = list(adapt_delta = 0.95),
file = "fits/time2",
file_refit = "on_change",
seed = 20210421
)
summary(time2)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: time ~ 1 + high_debt_version + scenario + mo(education_level) + (1 | session)
## Data: as.data.frame(d.both_completed) (Number of observations: 44)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~session (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.32 0.17 0.01 0.63 1.00 673 975
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 8.25 0.37 7.60 9.03 1.01 1415
## high_debt_versionfalse -0.15 0.16 -0.47 0.17 1.00 4332
## scenariotickets -0.32 0.17 -0.64 0.01 1.00 5152
## moeducation_level -0.24 0.12 -0.47 -0.01 1.01 1299
## Tail_ESS
## Intercept 1776
## high_debt_versionfalse 2746
## scenariotickets 2521
## moeducation_level 1913
##
## Simplex Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## moeducation_level1[1] 0.34 0.16 0.06 0.66 1.00 2900
## moeducation_level1[2] 0.16 0.11 0.02 0.45 1.00 3883
## moeducation_level1[3] 0.17 0.12 0.02 0.46 1.00 4003
## moeducation_level1[4] 0.32 0.15 0.07 0.63 1.00 4069
## Tail_ESS
## moeducation_level1[1] 2878
## moeducation_level1[2] 2719
## moeducation_level1[3] 2963
## moeducation_level1[4] 3137
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 3.97 1.24 2.11 6.91 1.00 947 2244
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
ranef(time2)
## $session
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## 6033d69a5af2c702367b3a95 -0.06489116 0.2373079 -0.56041770 0.4161218
## 6033d90a5af2c702367b3a96 0.14098737 0.2514511 -0.29310287 0.7132378
## 6034fc165af2c702367b3a98 0.24017793 0.2681683 -0.17918125 0.8398987
## 603500725af2c702367b3a99 -0.45162733 0.4026797 -1.31337375 0.1285293
## 603f97625af2c702367b3a9d -0.10064558 0.2423361 -0.62516282 0.3661429
## 603fd5d95af2c702367b3a9e 0.05787041 0.2444601 -0.42858840 0.5860447
## 60409b7b5af2c702367b3a9f 0.34139368 0.2947144 -0.08852588 0.9886314
## 604b82b5a7718fbed181b336 -0.13188201 0.2644422 -0.71980265 0.3506500
## 6050c1bf856f36729d2e5218 0.15362917 0.2503672 -0.27582955 0.7191141
## 6050e1e7856f36729d2e5219 -0.03175075 0.2450888 -0.54740960 0.4872661
## 6055fdc6856f36729d2e521b -0.08358514 0.2496640 -0.60834568 0.4079105
## 60589862856f36729d2e521f 0.14534943 0.2611730 -0.32717898 0.7414337
## 605afa3a856f36729d2e5222 -0.04418395 0.2387607 -0.53486885 0.4489215
## 605c8bc6856f36729d2e5223 -0.16713611 0.2752103 -0.79219780 0.3118615
## 605f3f2d856f36729d2e5224 -0.03562485 0.2365540 -0.51963285 0.4530460
## 605f46c3856f36729d2e5225 -0.06168182 0.2582907 -0.61540618 0.4564646
## 60605337856f36729d2e5226 0.21400610 0.2712890 -0.21215202 0.8405765
## 60609ae6856f36729d2e5228 0.21802292 0.2638524 -0.19673462 0.7920115
## 6061ce91856f36729d2e522e -0.12484399 0.2626830 -0.69784477 0.3467587
## 6061f106856f36729d2e5231 -0.33444313 0.3365499 -1.05775775 0.1775895
## 6068ea9f856f36729d2e523e 0.31627952 0.3397728 -0.20239922 1.0827567
## 6075ab05856f36729d2e5247 -0.11935109 0.2655445 -0.71282735 0.3687091
plot(time2, ask = FALSE)
pp_check(time2, nsamples = 200) + scale_x_log10()
We select the second best performing model with two variables.
time3 <- brm(
"time ~ 1 + high_debt_version + scenario + education_field + (1 | session)",
prior = c(
prior(normal(0, 1), class = "b"),
prior(normal(7.5, 1), class = "Intercept"),
prior(exponential(1), class = "sd"),
prior(gamma(0.01, 0.01), class = "shape")
),
family = negbinomial(),
data = as.data.frame(d.both_completed),
control = list(adapt_delta = 0.95),
file = "fits/time3",
file_refit = "on_change",
seed = 20210421
)
summary(time3)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: time ~ 1 + high_debt_version + scenario + education_field + (1 | session)
## Data: as.data.frame(d.both_completed) (Number of observations: 44)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~session (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.35 0.16 0.03 0.65 1.00 734 848
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept 7.65 0.23 7.21 8.11 1.00
## high_debt_versionfalse -0.17 0.16 -0.48 0.15 1.00
## scenariotickets -0.33 0.16 -0.65 -0.01 1.00
## education_fieldInteractionDesign -0.46 0.52 -1.47 0.57 1.00
## education_fieldNone 1.06 0.49 0.13 2.05 1.00
## education_fieldSoftwareEngineering 0.00 0.25 -0.50 0.47 1.00
## Bulk_ESS Tail_ESS
## Intercept 2651 2411
## high_debt_versionfalse 5421 2421
## scenariotickets 4774 2579
## education_fieldInteractionDesign 2879 2518
## education_fieldNone 2997 2649
## education_fieldSoftwareEngineering 2214 1912
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 4.14 1.27 2.19 6.95 1.00 1164 2395
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
ranef(time3)
## $session
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## 6033d69a5af2c702367b3a95 -0.08293346 0.2762997 -0.6506655 0.4595537
## 6033d90a5af2c702367b3a96 0.28430349 0.2753035 -0.1679328 0.8754226
## 6034fc165af2c702367b3a98 0.25582499 0.2663495 -0.1811189 0.8351591
## 603500725af2c702367b3a99 -0.43005665 0.3568812 -1.1641978 0.1529283
## 603f97625af2c702367b3a9d -0.13018561 0.2717640 -0.7127580 0.3790718
## 603fd5d95af2c702367b3a9e 0.18674536 0.2593710 -0.2540311 0.7503900
## 60409b7b5af2c702367b3a9f 0.36812731 0.2939810 -0.0971086 0.9932092
## 604b82b5a7718fbed181b336 -0.28082771 0.2978583 -0.9004351 0.2459909
## 6050c1bf856f36729d2e5218 0.15968857 0.2621604 -0.3126920 0.7298558
## 6050e1e7856f36729d2e5219 0.08134850 0.2504322 -0.3900459 0.6282933
## 6055fdc6856f36729d2e521b -0.04555245 0.2633609 -0.5669126 0.4875088
## 60589862856f36729d2e521f 0.01574282 0.2465412 -0.4567513 0.5440192
## 605afa3a856f36729d2e5222 -0.06195952 0.2684858 -0.6160360 0.4804998
## 605c8bc6856f36729d2e5223 -0.30838067 0.3069765 -0.9366980 0.2177001
## 605f3f2d856f36729d2e5224 -0.05390679 0.2672687 -0.5984670 0.4632304
## 605f46c3856f36729d2e5225 -0.20080961 0.2769727 -0.8019196 0.3072222
## 60605337856f36729d2e5226 0.22553691 0.2748839 -0.2484014 0.8242951
## 60609ae6856f36729d2e5228 0.30829846 0.2922382 -0.1723765 0.9343437
## 6061ce91856f36729d2e522e -0.02434554 0.2449203 -0.5018119 0.4834777
## 6061f106856f36729d2e5231 -0.28477542 0.3054150 -0.9283804 0.2438914
## 6068ea9f856f36729d2e523e 0.12585163 0.3566877 -0.5549790 0.9069138
## 6075ab05856f36729d2e5247 -0.07493014 0.3733400 -0.8785455 0.6827192
plot(time3, ask = FALSE)
pp_check(time3, nsamples = 200) + scale_x_log10()
All candidate models look nice, none is significantly better than the others, we will proceed the simplest model: time0
We will try a few different variations of the selected candidate model.
Some participants did only complete one scenario. Those has been excluded from the initial dataset to improve sampling of the models. We do however want to use all data we can and will therefore try to fit the model with the complete dataset.
time0.all <- brm(
"time ~ 1 + high_debt_version + (1 | session)",
prior = c(
prior(normal(0, 1), class = "b"),
prior(normal(7.5, 1), class = "Intercept"),
prior(exponential(1), class = "sd"),
prior(gamma(0.01, 0.01), class = "shape")
),
family = negbinomial(),
data = as.data.frame(d.completed),
control = list(adapt_delta = 0.95),
file = "fits/time0.all",
file_refit = "on_change",
seed = 20210421
)
summary(time0.all)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: time ~ 1 + high_debt_version + (1 | session)
## Data: as.data.frame(d.completed) (Number of observations: 51)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~session (Number of levels: 29)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.43 0.14 0.13 0.69 1.00 926 958
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 7.53 0.14 7.26 7.82 1.00 2362
## high_debt_versionfalse -0.20 0.16 -0.49 0.11 1.00 4301
## Tail_ESS
## Intercept 2409
## high_debt_versionfalse 2972
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 3.76 1.05 2.08 6.08 1.00 1169 1896
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
ranef(time0.all)
## $session
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## 6033c6fc5af2c702367b3a93 0.10743406 0.3215793 -0.49176255 0.7797215
## 6033d69a5af2c702367b3a95 -0.11393560 0.2878708 -0.68202983 0.4603690
## 6033d90a5af2c702367b3a96 0.32426832 0.2884523 -0.18827020 0.9204864
## 6034fc165af2c702367b3a98 0.32566422 0.2880610 -0.19423385 0.9328176
## 603500725af2c702367b3a99 -0.51453087 0.3640134 -1.24499075 0.1281945
## 603f84f15af2c702367b3a9b -0.36670177 0.3991751 -1.18778925 0.3558798
## 603f97625af2c702367b3a9d -0.18153951 0.2919422 -0.74118115 0.4041113
## 603fd5d95af2c702367b3a9e 0.25168127 0.2876644 -0.27646243 0.8557478
## 60409b7b5af2c702367b3a9f 0.38228571 0.3033029 -0.16228467 1.0108630
## 604b82b5a7718fbed181b336 -0.33134416 0.3220736 -0.98139030 0.2840337
## 604f1239a7718fbed181b33f -0.02585933 0.3359783 -0.66638825 0.6747302
## 6050c1bf856f36729d2e5218 0.12702070 0.2844274 -0.39857212 0.7282526
## 6050e1e7856f36729d2e5219 0.08732978 0.2788988 -0.42141037 0.6620035
## 6055fdc6856f36729d2e521b -0.06830187 0.2847831 -0.63595455 0.5284143
## 60579f2a856f36729d2e521e -0.05001702 0.3367932 -0.70177923 0.6469537
## 60589862856f36729d2e521f 0.05826002 0.2851661 -0.47828292 0.6554365
## 605a30a7856f36729d2e5221 -0.20497733 0.3589553 -0.95256537 0.4872980
## 605afa3a856f36729d2e5222 -0.06369539 0.2935347 -0.64654803 0.5295743
## 605c8bc6856f36729d2e5223 -0.36294381 0.3313016 -1.02561925 0.2553571
## 605f3f2d856f36729d2e5224 -0.11411407 0.2794536 -0.64467565 0.4437461
## 605f46c3856f36729d2e5225 -0.25334110 0.3045431 -0.85241390 0.3196846
## 60605337856f36729d2e5226 0.24639023 0.2910492 -0.26199735 0.8719473
## 60609ae6856f36729d2e5228 0.33274912 0.2890039 -0.18430892 0.9341700
## 6061ce91856f36729d2e522e -0.03236292 0.2815524 -0.57206517 0.5407210
## 6061f106856f36729d2e5231 -0.37393379 0.3343367 -1.00906150 0.2773342
## 60672faa856f36729d2e523c -0.12075577 0.3383869 -0.79444963 0.5510030
## 6068ea9f856f36729d2e523e 0.69985604 0.3314307 0.03456946 1.3554313
## 606db69d856f36729d2e5243 0.57103782 0.3699743 -0.05728208 1.3535968
## 6075ab05856f36729d2e5247 -0.30871951 0.3183636 -0.92231982 0.2797401
plot(time0.all, ask = FALSE)
pp_check(time0.all, nsamples = 200) + scale_x_log10()
As including all data points didn’t harm the model we will create this variant with all data points as well.
time0.all.exp <- brm(
"time ~ 1 + high_debt_version + work_experience_programming.s + (1 | session)",
prior = c(
prior(normal(0, 1), class = "b"),
prior(normal(7.5, 1), class = "Intercept"),
prior(exponential(1), class = "sd"),
prior(gamma(0.01, 0.01), class = "shape")
),
family = negbinomial(),
data = as.data.frame(d.completed),
control = list(adapt_delta = 0.95),
file = "fits/time0.all.exp",
file_refit = "on_change",
seed = 20210421
)
summary(time0.all.exp)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: time ~ 1 + high_debt_version + work_experience_programming.s + (1 | session)
## Data: as.data.frame(d.completed) (Number of observations: 51)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~session (Number of levels: 29)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.45 0.14 0.14 0.71 1.01 708 654
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept 7.54 0.14 7.26 7.82 1.00
## high_debt_versionfalse -0.20 0.16 -0.50 0.10 1.00
## work_experience_programming.s 0.04 0.12 -0.20 0.28 1.00
## Bulk_ESS Tail_ESS
## Intercept 1883 2616
## high_debt_versionfalse 4496 2962
## work_experience_programming.s 1589 2248
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 3.75 1.05 2.04 6.02 1.00 1059 1646
##
## Samples were drawn using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
ranef(time0.all.exp)
## $session
## , , Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## 6033c6fc5af2c702367b3a93 0.11436634 0.3386521 -0.52261055 0.8266039
## 6033d69a5af2c702367b3a95 -0.10323867 0.2966430 -0.68971485 0.4861681
## 6033d90a5af2c702367b3a96 0.33873191 0.2935511 -0.18501772 0.9456097
## 6034fc165af2c702367b3a98 0.34495015 0.3037295 -0.19167927 0.9722129
## 603500725af2c702367b3a99 -0.51896554 0.3736859 -1.26467925 0.1508850
## 603f84f15af2c702367b3a9b -0.38313522 0.4158486 -1.23419975 0.3778444
## 603f97625af2c702367b3a9d -0.17090340 0.3091673 -0.76201672 0.4497196
## 603fd5d95af2c702367b3a9e 0.27638045 0.2904963 -0.26957063 0.8767928
## 60409b7b5af2c702367b3a9f 0.40718387 0.3049947 -0.13132455 1.0372023
## 604b82b5a7718fbed181b336 -0.33476385 0.3332111 -0.98945650 0.2977670
## 604f1239a7718fbed181b33f -0.03143005 0.3278864 -0.67265365 0.6432929
## 6050c1bf856f36729d2e5218 0.12539732 0.2860831 -0.42352197 0.7080399
## 6050e1e7856f36729d2e5219 0.08576607 0.2765233 -0.43495113 0.6692808
## 6055fdc6856f36729d2e521b -0.05793969 0.2966705 -0.64860180 0.5421729
## 60579f2a856f36729d2e521e -0.05592628 0.3547309 -0.75353255 0.6811279
## 60589862856f36729d2e521f 0.01215595 0.3402376 -0.65264142 0.7036291
## 605a30a7856f36729d2e5221 -0.20484817 0.3679164 -0.96836420 0.5002716
## 605afa3a856f36729d2e5222 -0.10318751 0.3136994 -0.72476087 0.5194962
## 605c8bc6856f36729d2e5223 -0.38926449 0.3301529 -1.05361050 0.2138376
## 605f3f2d856f36729d2e5224 -0.17747834 0.3823643 -0.98431508 0.5444238
## 605f46c3856f36729d2e5225 -0.24906402 0.3145777 -0.89621808 0.3507401
## 60605337856f36729d2e5226 0.26867787 0.2941246 -0.24631512 0.8860694
## 60609ae6856f36729d2e5228 0.36009662 0.2998216 -0.17088095 0.9726449
## 6061ce91856f36729d2e522e -0.02747183 0.2857916 -0.57670435 0.5517471
## 6061f106856f36729d2e5231 -0.37962825 0.3406161 -1.04951750 0.2722603
## 60672faa856f36729d2e523c -0.11292912 0.3609199 -0.82503188 0.6248903
## 6068ea9f856f36729d2e523e 0.71307170 0.3365822 0.03128014 1.3653205
## 606db69d856f36729d2e5243 0.55498432 0.3645415 -0.07069035 1.3209760
## 6075ab05856f36729d2e5247 -0.30613297 0.3229255 -0.95367582 0.2899334
loo(
time0.all,
time0.all.exp
)
## Output of model 'time0.all':
##
## Computed from 4000 by 51 log-likelihood matrix
##
## Estimate SE
## elpd_loo -425.1 7.8
## p_loo 17.5 3.1
## looic 850.2 15.6
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 25 49.0% 758
## (0.5, 0.7] (ok) 21 41.2% 258
## (0.7, 1] (bad) 5 9.8% 70
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'time0.all.exp':
##
## Computed from 4000 by 51 log-likelihood matrix
##
## Estimate SE
## elpd_loo -425.7 7.8
## p_loo 18.2 3.2
## looic 851.4 15.5
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 30 58.8% 756
## (0.5, 0.7] (ok) 16 31.4% 362
## (0.7, 1] (bad) 4 7.8% 42
## (1, Inf) (very bad) 1 2.0% 12
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## time0.all 0.0 0.0
## time0.all.exp -0.6 0.5
plot(time0.all.exp, ask = FALSE)
pp_check(time0.all.exp, nsamples = 200) + scale_x_log10()
This means that our final model, with all data points and experience predictors, is time0.all.exp
To begin interpreting the model we look at how it’s parameters were estimated. As our research is focused on how the outcome of the model is effected we will mainly analyze the \(\beta\) parameters.
mcmc_areas(time0.all.exp, pars = c("b_high_debt_versionfalse", "b_work_experience_programming.s"), prob = 0.95) + scale_y_discrete() +
scale_y_discrete(labels=c("High debt version: false", "Professional programming experience")) +
ggtitle("Beta parameters densities in time model", subtitle = "Shaded region marks 95% of the density. Line marks the median")
We start by extracting posterior samples
scale_programming_experience <- function(x) {
(x - mean(d.completed$work_experience_programming))/ sd(d.completed$work_experience_programming)
}
unscale_programming_experience <- function(x) {
x * sd(d.completed$work_experience_programming) + mean(d.completed$work_experience_programming)
}
post_settings <- expand.grid(
high_debt_version = c("false", "true"),
session = NA,
work_experience_programming.s = sapply(c(0, 3, 10, 25, 40), scale_programming_experience)
)
post <- posterior_predict(time0.all.exp, newdata = post_settings) %>%
melt(value.name = "estimate", varnames = c("sample_number", "settings_id")) %>%
left_join(
rowid_to_column(post_settings, var= "settings_id"),
by = "settings_id"
) %>%
mutate(work_experience_programming = unscale_programming_experience(work_experience_programming.s)) %>%
select(
estimate,
high_debt_version,
work_experience_programming
) %>%
mutate(estimate = estimate/60)
ggplot(post, aes(x=estimate, fill = high_debt_version)) +
geom_density(alpha = 0.5) +
scale_x_log10() +
scale_fill_manual(
name = "Debt version",
labels = c("Low debt", "High debt"),
values = c("lightblue", "darkblue")
) +
facet_grid(rows = vars(work_experience_programming)) +
labs(
title = "Time to complete task / years of programming experience",
subtitle = "Notice! x-axis is log10 scaled.",
x = "Time (min)",
y = "Density"
)
post.diff <- post %>% filter(high_debt_version == "true")
post.diff$estimate = post.diff$estimate - filter(post, high_debt_version == "false")$estimate
ggplot(post.diff, aes(x=estimate, y = 0, fill = stat(quantile))) +
geom_density_ridges_gradient(quantile_lines = TRUE, quantile_fun = hdi, vline_linetype = 2) +
scale_fill_manual(name = "HDI", labels = c("100%", "95%", "100%"), values = c("transparent", "lightblue", "transparent"),) +
xlim(-100, 100) +
facet_grid(rows = vars(work_experience_programming)) +
labs(
title = "Time to complete task / years of programming experience",
subtitle = "Difference as: high debt time - low debt time",
x = "Time (min)",
y = "Density"
)